* By Hsin-Wei Su and Hung-Jen Wang

* version 1.0  2012/09/02

capture program drop sfsicm
program define sfsicm

version 10.1

syntax varlist, Distribution(string) FRONTIER(string) Brep(string) /*
                */ [COST PRODuction  Stat(string) C(real 5) MU(string) USIGMAS(string) VSIGMAS(string) /*
                */ ETAS(string) INIT(string) SEED(string) NODOTS]

local themat = `brep'+ 10
quie sum `varlist'
local totalN = r(N)  /* to be used later */
local thematb = r(N) + 10

local largemat = max(`themat', `thematb')

set matsize `largemat'

capture set seed `seed'

unab frontier: `frontier'
capture unab mu: `mu'
capture unab usigmas: `usigmas'
capture unab vsigmas: `vsigmas'
capture unab etas: `etas'

*---- Check if the exogenous vector is the same -----------*

local wrong = 0
if ("`mu'"~="") & ("`usigmas'" ~= "") & ("`mu'" ~= "`usigmas'") {
   local wrong = 1
}
if ("`mu'"~="") & ("`vsigmas'" ~= "") & ( "`mu'" ~= "`vsigmas'") {
   local wrong = 1
}
if ("`usigmas'"~="") & ("`vsigmas'" ~= "") & ( "`usigmas'" ~= "`vsigmas'") {
   local wrong = 1
}
if ("`mu'"~="") & ("`etas'" ~= "")  {
   di in red "You can't have both -mu- (for truncated normal) and -etas- (for exponential) specified together."
   error 199
   exit
}
if ("`etas'"~="") & ("`vsigmas'" ~= "") & ( "`etas'" ~= "`vsigmas'") {
   local wrong = 1
}
if `wrong' == 1 {
 di in red "The exogenous variables in mu, usigmas, and vsigmas need to be the same vector."
 error 199
 exit
}

*------ setup zlist and initial values ------

if ("`mu'"~="") {
  local zlist `mu'
}
else if ("`usigmas'" ~= "") {
  local zlist `usigmas'
}
else if ("`vsigmas'" ~= "") {
  local zlist `vsigmas'
}
else if ("`etas'" ~= "") {
  local zlist `etas'
}

matrix im=0
foreach M of local mu {
     matrix im=im\0
}
matrix iu=0
foreach U of local usigmas {
     matrix iu=iu\0
}
matrix iv=0
foreach V of local vsigmas {
     matrix iv=iv\0
}
matrix ie=0
foreach E of local etas {
     matrix ie=ie\0
}

* -------------------------


tempvar oldid
gen `oldid' = _n

*-------------------------------------
global PorC = "Correct PorC Specification" /* offset previous values, so that errors are issued if a proper value is not given */

if "`production'" ~= "" & "`cost'" ~= "" {
   di in red "You can specify only one of the -cost- or -prod-, not both."
}

if "`production'" ~= "" {
   global PorC = 1
}
else if "`cost'" ~= ""{
   global PorC = 2
}
else {
   di in red "You need to specify -cost- or -production-."
   exit 198
}
* ------------------------------------

local needinit = 1
if "`init'" ~= "" { /* user supplied the initial values */
 local needinit = 0
}


capture reg `varlist' `frontier'
mat b0 = e(b)

if  ("`distribution'" == "h") | ("`distribution'" == "halfnormal") {
  sfmodel `varlist', `production' `cost' dist(h) frontier(`frontier') usigmas(`usigmas') vsigmas(`vsigmas')
  if `needinit' == 1 {
    sf_init, frontier(b0) usigmas(iu) vsigmas(iv)
  }
}

if ("`distribution'" == "t") | ("`distribution'" == "truncated") {
  sfmodel `varlist', `production' `cost' dist(t) frontier(`frontier') mu(`mu') usigmas(`usigmas') vsigmas(`vsigmas')
  if `needinit' == 1 {
     sf_init, frontier(b0) mu(im) usigmas(iu) vsigmas(iv)
  }
}
if ("`distribution'" == "e") | ("`distribution'" == "exponential") {
  sfmodel `varlist', `production' `cost' dist(e) frontier(`frontier') etas(`etas') vsigmas(`vsigmas')
  if `needinit' == 1 {
    sf_init, frontier(b0) etas(ie) vsigmas(iv)
  }
}

  if `needinit' == 1 {
    sf_srch, n(2) frontier(`frontier') mu(`mu') usigmas(`usigmas') vsigmas(`vsigmas') etas(`etas') fast nograph
  }
  else {
   ml init `init', copy
  }

   ml max, diff gtol(1e-4) nrtol(1e-4) iterate(50)

scalar nofit = e(ic)

global badrep = 0 /* This line is mainly for simulation purpose. */

if  e(ic) == 50 { /* meaning the model does not converge smoothly */
  di " "
  di in red "The sf model does not converge smoothly."
  di " "
  global badrep = 1 /* This line is mainly for simulation purpose. */
  error 1999
}

matrix mb=e(b)'

tempvar ty tu tv pf pmu pusigln pvsigln us vs es rdn petasln

    if ("`distribution'" == "e") | ("`distribution'" == "exponential") {
        predict double `petasln', eq(etas) xb
        gen double `es' = exp(`petasln')
        gen double `rdn' = uniform()
        gen double `tu' = -(sqrt(`es'))*ln(1-`rdn')
    }
    else { /* half normal or truncated normal */
     if ("`distribution'" == "t") | ("`distribution'" == "truncated") {
        predict double `pmu', eq(mu) xb
     }
     if ("`distribution'" == "h") | ("`distribution'" == "halfnormal") {
        gen `pmu' = 0
     }
      predict double `pusigln', eq(usigmas) xb
      gen double `us' = exp(`pusigln')
      gentrun double `tu', dist(`pmu' `us') left(0)
      quie replace `tu' = 0 if `tu' == .
    }


 predict double `pvsigln', eq(vsigmas) xb
 gen double `vs' = exp(`pvsigln')
 gentrun double `tv', dist(0 `vs')

 predict double `pf', eq(frontier) xb


if $PorC == 2 {
  quie replace `tu' = -`tu'
}

gen double `ty'=`pf'+`tv'-`tu'


tempvar ys yt ay ayt
egen double `ys'=std(`varlist')
egen double `yt'=std(`ty')
gen double `ay'=atan(`ys')
gen double `ayt'=atan(`yt')

drop `tu' `tv' `ty' `ys' `yt'


*foreach X of local zlist {  /* to avoid overlpping varaibles in xvar and zlist, change the name of zlist */
*  tempvar z`X'
*  quie gen double `z`X'' = `X'
*  local zzlist "`z`X''"
*}

*local zlist `zzlist'

local xzvar = "`frontier' `zlist'"
local xzvar: list uniq xzvar

foreach X of local xzvar {
     tempvar std_`X'
     egen double `std_`X''=std(`X')
}

foreach X of local xzvar {
     tempvar atan_`X'
     gen double `atan_`X''=atan(`std_`X'')
     drop `std_`X''
}

local N=_N
local N_1=_N-1

tempname eps s0 s1 s2
scalar `eps' =10^(-100)
* scalar c=`c'
scalar `s1'=0

tempvar px ay_i ayt_i bf1 bf2

foreach X of local xzvar { /* create the temp var for later use */
  tempvar ax_i`X'
}

forvalues i=1(1)`N_1' {
  tempvar mpx`i'
}

quietly {
forvalues i=1(1)`N_1' {

     gen double `ay_i'=`ay'[`i']
     gen double `ayt_i'=`ayt'[`i']
     local i_1=`i'+1
     gen double `px'=1
     foreach X of local xzvar {
          gen double `ax_i`X''=`atan_`X''[`i']
     }

     foreach X of local xzvar {
          gen double `bf1' = `px'
          drop `px'
          gen double `px'=`bf1'*sin((`c')*(`ax_i`X''-`atan_`X''+scalar(`eps')))/((`c')*(`ax_i`X''-`atan_`X''+scalar(`eps')))
          drop `ax_i`X'' `bf1'
     }

     egen double `bf2' = total((sin((`c')*(`ay_i'-`ay'+scalar(`eps')))/((`c')*(`ay_i'-`ay'+scalar(`eps')))-sin((`c')*(`ayt_i'-`ay'+scalar(`eps')))/((`c')*(`ayt_i'-`ay'+scalar(`eps')))-sin((`c')*(`ay_i'-`ayt'+scalar(`eps')))/((`c')*(`ay_i'-`ayt'+scalar(`eps')))+sin((`c')*(`ayt_i'-`ayt'+scalar(`eps')))/((`c')*(`ayt_i'-`ayt'+scalar(`eps'))))*`px') in `i_1'/`N'

     scalar `s0'=`bf2'[`N']

     capture drop `mpx`i''
     gen double `mpx`i''=`px'

     drop `ay_i' `ayt_i' `px' `bf2'
     scalar `s1' =scalar(`s1')+scalar(`s0')
}

tempvar bf3
gen double `bf3'=sum(sin((`c')*(`ay'-`ayt'+scalar(`eps')))/((`c')*(`ay'-`ayt'+scalar(`eps'))))
* tempname s2
scalar `s2'=`bf3'[`N']

scalar sicm_test_stat=2+2/`N'*(scalar(`s1')-scalar(`s2'))
drop `ay' `ayt' `bf3'
}

if "`stat'" ~="" {
  scalar `stat'=sicm_test_stat
}



di " "
di in yel "Null: The model specification is correct."
di in yel "The SICM test statistic is " sicm_test_stat "."
di in yel "     Now bootstrapping the distribution... Be patient."
di " "

matrix SICM_dist=0

if "`nodots'" == "" {
  _dots 0, title(bootstrap running) reps(`brep')
}

tempvar sy su sv rdnn
tempvar ty tu tv tpf tmz bpusigln uss bpvsigln vss bpetasln bes rdnnn
tempname s0 s2


forvalues p=1(1)`brep' {

*     tempvar sy su sv rdnn
     matrix mbs=mb

    if ("`distribution'" == "e") | ("`distribution'" == "exponential") {
        gen double `rdnn' = uniform()
        gen double `su' = -(sqrt(`es'))*ln(1-`rdnn')
        drop `rdnn'
    }
    else { /* half or truncated normal */
        gentrun double `su', dist(`pmu' `us') left(0)
        quie replace `su' = 0 if `su' == .
    }

     gentrun double `sv', dist(0 `vs')

     if $PorC == 2 {
       quie replace `su' = -`su'
     }

     gen double `sy'=`pf'+`sv'-`su'

     drop `su' `sv'


       if  ("`distribution'" == "h") | ("`distribution'" == "halfnormal") {
         sfmodel `sy', `production' `cost' dist(h) frontier(`frontier') usigmas(`usigmas') vsigmas(`vsigmas')
       }

       if ("`distribution'" == "t") | ("`distribution'" == "truncated") {
         sfmodel `sy', `production' `cost' dist(t) frontier(`frontier') mu(`mu') usigmas(`usigmas') vsigmas(`vsigmas')
       }
      if ("`distribution'" == "e") | ("`distribution'" == "exponential") {
         sfmodel `sy', `production' `cost' dist(e) frontier(`frontier') etas(`etas') vsigmas(`vsigmas')
       }

         ml init mb, copy
         capture ml max, difficult nonrtolerance tol(1e-4) iterate(50) search(off)

     if _rc == 0 & e(ic) < 50 {

     matrix mbs=e(b)'

*     tempvar ty tu tv tpf tmz bpusigln uss bpvsigln vss bpetasln bes rdnnn

    if ("`distribution'" == "e") | ("`distribution'" == "exponential") {
        predict double `bpetasln', eq(etas) xb
        gen double `bes' = exp(`bpetasln')
        gen double `rdnnn' = uniform()
        gen double `tu' = -(sqrt(`bes'))*ln(1-`rdnnn')
    }
    else { /* half normal or truncated normal */
      if ("`distribution'" == "t") | ("`distribution'" == "truncated") {
        predict double `tmz', eq(mu) xb
      }
      if ("`distribution'" == "h") | ("`distribution'" == "halfnormal") {
        gen `tmz' = 0
      }

     predict double `bpusigln', eq(usigmas) xb
     gen double `uss' = exp(`bpusigln')
     gentrun double `tu', dist(`tmz' `uss') left(0)
     quie replace `tu' = 0 if `tu' == .

    }

     predict double `bpvsigln', eq(vsigmas) xb
     gen double `vss' = exp(`bpvsigln')
     gentrun double `tv', dist(0 `vss')

     predict double `tpf', eq(frontier) xb

     if $PorC == 2 {
       quie replace `tu' = -`tu'
     }

     gen double `ty'=`tpf'+`tv'-`tu'

     egen double `ys'=std(`sy')
     egen double `yt'=std(`ty')
     gen double `ay'=atan(`ys')
     gen double `ayt'=atan(`yt')

     drop `tu' `tv' `tpf' `vss'
     capture drop `tmz'
     capture drop `uss'
     capture drop `bpetasln'
     capture drop `bes'
     capture drop `rdnnn'

     capture drop `bpusigln'
     capture drop `bpvsigln'

     drop `sy' `ty' `ys' `yt'

     scalar `s1'=0

     quietly {
     forvalues i=1(1)`N_1' {
          gen double `ay_i'=`ay'[`i']
          gen double `ayt_i'=`ayt'[`i']
          local i_1=`i'+1
          egen double `bf2' = total((sin((`c')*(`ay_i'-`ay'+scalar(`eps')))/((`c')*(`ay_i'-`ay'+scalar(`eps')))-sin((`c')*(`ayt_i'-`ay'+scalar(`eps')))/((`c')*(`ayt_i'-`ay'+scalar(`eps')))-sin((`c')*(`ay_i'-`ayt'+scalar(`eps')))/((`c')*(`ay_i'-`ayt'+scalar(`eps')))+sin((`c')*(`ayt_i'-`ayt'+scalar(`eps')))/((`c')*(`ayt_i'-`ayt'+scalar(`eps'))))*`mpx`i'') in `i_1'/`N'
*          tempname s0
          scalar `s0'=`bf2'[`N']
          drop `ay_i' `ayt_i' `bf2'

          scalar `s1'=scalar(`s1')+scalar(`s0')
     }
     gen double `bf3'=sum(sin((`c')*(`ay'-`ayt'+scalar(`eps')))/((`c')*(`ay'-`ayt'+scalar(`eps'))))
*     tempname s2
     scalar `s2'=`bf3'[`N']
     scalar sicm_test_stat=2+2/`N'*(scalar(`s1')-scalar(`s2'))

     drop `ay' `ayt' `bf3'
     }
     matrix SICM_dist=SICM_dist\sicm_test_stat


     }  /* if _rc == 0 */
   else { /* if the replication does not converge */

     foreach X in `tu' `tv' `tpf' `vss' `tmz' `uss' `bpetasln' `bes' `rdnnn' `bpusigln' `bpvsigln' {
        capture drop `X'
     }

    foreach X in `sy' `ty' `ys' `yt' `ay_i' `ayt_i' `bf2' `ay' `ayt' `bf3' {
        capture drop `X'
    }

   }

    if "`nodots'"  == "" {
        _dots `p' 0
     }

*    scalar qqq = 10*floor(`p'/10)
*    if qqq == `p' {
*       memory
*     }

} /* bootstrap */


matrix SICM_dist=SICM_dist[2...,1]

capture drop SICM_dist1
svmat SICM_dist

quie summarize SICM_dist1, detail

tempname r90 r95 r99

scalar `r90' = r(p90)
scalar `r95' = r(p95)
scalar `r99' = r(p99)

/* The following 4 lines are mainly for simulation purpose. */
scalar r90 = `r90'
scalar r95 = `r95'
scalar r99 = `r99'
scalar nofboot = rowsof(SICM_dist)



drop SICM_dist1
sort `oldid'
keep in 1/`totalN'

di " "
di in yel "Values of the bootstrapped distribution is store in the matrix SICM_dist."
di " "
di in yel "The 1%, 5%, and 10% quantile values are, respectively, "
di in yel "     " `r99' ", "
di in yel "     " `r95' ", and "
di in yel "     " `r90' "."


end
